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We propose a formalism to study dynamical properties of a quantum many-body system in the 
thermodynamic limit by studying a finite system with "infinite boundary conditions" (IBC) where 
both finite size effects and boundary effects have been eliminated. For one-dimensional systems, 
infinite boundary conditions are obtained by attaching two boundary sites to a finite system, where 
each of these two sites effectively represents a semi-infinite extension of the system. One can then 
use standard finite-size matrix product state techniques to study a region of the system while avoid- 
ing many of the complications normally associated with finite-size calculations such as boundary 
Friedel oscillations. We illustrate the technique with an example of time evolution of a local per- 
turbation applied to an infinite (translationally invariant) ground state, and use this to calculate 
the spectral function of the 5=1 Heisenberg spin chain. This approach is more efficient and more 
accurate than conventional simulations based on finite-size matrix product state and density-matrix 
renormalization-group approaches. 
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I. INTRODUCTION 

In recent decades, the tensor network formalism has 
emerged as a set of powerful numerical techniques to in- 
vestigate physical properties of strongly-correlated quan- 
tum many-body systems. For instance, in ID sys tem s, 
the density-matrix renormalization-group (DMRG^|i^ is 
probably the single most powerful method to compute 
numerically exact ground states. Furthermore, the de- 
velopment of the time-evolving block decimation (TEBD) 



algorithm'^ highlighted the great advantages of the ma- 
trix product state (MPS^^P^^ representation, which incor- 
porates DMRG and TEBD into the same frameworlP^l. 
Meanwhile, tensor product state (TPSJP^HUl and pro- 
jected entangled-pair state (PEPSJiSHiIl methods are de- 
veloping into important tools for the study of 2D systems. 

For calculating bulk properties of matter, it is desirable 
to take the thermodynamic limit and avoid the influence 
of boundary conditions. In many methods, the thermo- 
dynamic limit is not possible to study directly, but in- 
stead requires the extrapolation of results for increasingly 
larger system sizes. This is because for most algorithms 
the computational cost increases with the system size, 
however approaching the thermodynamic limit in this 
way is computationally expensive. In ID, there exist al- 
gorithms that overcome this limit by taking advantage of 
the invariance under translation in space. One of t hese is 
the infinite time-evolving block decimation (iTEBD) ^^ l ^° l, 
originally introduced to investigate the time evolution 
problem for infinite-size ID spin chains. In this algo- 
rithm, the infinite MPS (IMPS) is represented by a small 
set of tensors which are invariant under translation of one 
unit cell (equal to two sites for the usual TEBD scheme) . 
This algorithm can be used to obtain a translationally 
invariant ground state by evolving the tensors in imagi- 
nary time until the fixed point is reached. The resulting 
IMPS is not only a good representation of ground state, 
but compared with finite MPS the number of wavefunc- 
tion parameters is reduced and the IMPS form is very 
convenient for calculating observables of the system in 
the thermodynamic limit. The iTEBD algorithm is very 
easy to implement, however there are many ways to op- 
timize an IMPS to achieve the same fixed point. A faster 
converging algorithm which also allows more flexibility 
in the size of the unit cell is the iDMRG^^ algorithm, 
but other algorithms exist with some advantages for some 
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situationflSISl. 

Although the iMPS representation of a wavefunction 
is very useful for studying physical systems in the ther- 
modymanic limit, there are some applications for which 
breaking of translational invariance is essential, such as 
the response to a local perturbation. The time evolu- 
tion of a local perturbation is a common technique used 
in MPS calculations to obtain the spectral functiorP^H 
which to date has required using a finite MPS represen- 
tation. However, the use of a finite MPS has several 
disadvantages. In particular, the system size needs to 
be large enough that the excitation is not influenced by 
the boundary of the system. This clearly requires that 
the propogating excitation will not hit the boundary, but 
even this is not enough since the boundary will induce 
inhomogeneities such as Friedel oscillations, which means 
that the system size must be quite large even to obtain an 
approximately homogeneous ground state in the central 
region of the lattice. 

The notion of translational invariance of an IMPS can 
be generalized to states with finite momentum, whereby 
instead of requiring invariance under some number of lat- 
tice shifts, we instead require only that the IMPS is an 
eigenstate of translations with some complex eigenvalue 
gzfe representing the momentum. The resulting IMPS re- 
mains position independent but is constructed in such 
a way that the transfer operator has non-trivial phase 
factors. Algo rithms have been proposed for expecta- 
tion value^S^ and quasi-particle excitation^^Zl using this 
scheme. 

For infinite-size systems an equivalent problem was 
also investigated in Ref. [251 In that work the authors 
proposed an efficient method to simulate both imaginary- 
and real-time evolution of the infinite-size system with 
impurities by transversely contracting the tensor net- 
works along the space direction rather than along the 
time direction as in standard iTEBD. By using a folding 
technique to reduce the entangment of the MPS represen- 
tation the transverse contraction approach can achieve 
longer times than other techniques, nevertheless it can- 
not avoid some drawbacks. For instance, by employing 
the Suzuki- Trotter decompositiorP^ in the evolution op- 
erator with small time step, the finite number of rows 
along the time axis may be very large. This may cause 
difficulty in finding the left and right dominant eigenvec- 
tors of the transfer matrix. 

We will investigate the above problem in a different 
way by introducing what we call infinite boundary condi- 
tions for a finite MPS. We begin with the ground state of 
a many-body ID system described by an iMPS. A finite 
region of the infinite system can be perturbed while still 
utilizing the IMPS structure for the tensors not directly 
affected by the perturbation. The resulting structure is 
equivalent to a finite MPS with a specially constructed 
'pseudosite' at each end which effectively represents an 
infinite extension of the system. A key point of this con- 
struction is that the Hilbert space for the infinite exten- 
sion is fixed but the wavefunction is not, hence it can 



freely explore all of the available states in the effective 
Hilbert space of the infinite extension. The result is that, 
in contrast to conventional finite size MPS calculations 
where a propagating excitation reaching the boundary of 
the system will reflect back, in our scheme an excitation 
can propogate off the end of the finite MPS. As long as 
the perturbation outside the finite boundary is not too 
big there is little loss in fidelity from allowing it to do so. 

The evolution of a finite section of an IMPS was con- 
sidered by Kjall et aP^, who used this notion in obtain- 
ing the time evolution of a translationally invariant state 
that had been perturbed by a local particle excitation. 
However their scheme was rather specific to the particu- 
lar setting, of Suzuki-Trotter-based real-time evolution. 
In this paper we show that this idea can be taken much 
further, and by mapping the problem onto a finite MPS 
then any algorithm for finite MPS calculations can be 
applied to an infinite system. 

The paper is organised as follows: In Sec. II we will in- 
troduce the infinite boundary condition definition and ef- 
fective Hamiltonian calculation. In Sec. HI we review the 
problem of a local perturbation in the infinite spin chain 
and real-time evolution algorithm. In Sec. IV we then 
apply the idea of infinite boundary conditions to simulate 
the time evolution of spin-1 isotropic antiferromagnetic 
Heisenberg model. The results are presented by calcu- 
lating time-dependent observables such as local magne- 
tization {Sz{x,t)) to see how a wavefront propagates in 
time and unequal time two-point correlator A{x, t), from 
which we can extract the spectral function and disper- 
sion relation of the system. Finally, Sec. V contains our 
conclusions. 



II. INFINITE BOUNDARY CONDITIONS 

A. Formulation 

Let us consider an infinite-size spin chain for which the 
wave function is described by a one-site translationally 
invariant canonical iMPS 

I*) = ^ . . . Ar"'-iAr^'Ar"'+iAr'''+= . . . |s), (i) 

where |s) = | . . . Si-i, Si, Si+i, Si+2 . . .), Si is the local in- 
dex that represents an element in local Hilbert space at 
the ith site of the spin chain. The matrices V and A have 
dimension xx x and A is diagonal. Notice that while the 
notation of bond dimension is usually used as m or D 
in the DMRG language, here we use x instead, x plays 
a role as the refinement parameter of the IMPS. Specifi- 
cally, the larger the x the better the iMPS can represent 
the state. Diagrammatically, the iMPS is illustrated in 
Fig. [T^ where a pair of tensors {F, A} is repeated at every 
lattice site throughout the whole infinite chain. 

For convenience and later use, we can also rewrite Eq.[l] 
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FIG. 1. (Color online): (a). Single site translationally invari- 
ant canonical form IMPS. (b). One-site translationally in- 
variant mixed canonical form IMPS. (c). Partition the whole 
chain into three parts: the left and right semi-infinite sub- 
lattices, and the middle part which is a window that con- 
tains A'' sites, (d). Finite-size MPS effectively represent the 
IMPS with left and right-effective sites representing the left 
and right semi-infinite sublattices. 



represent the whole left and right semi-infinite chains, 
where each has dimension x ^ X? Fig. [l}i. These two 
matrices represent two boundary sites attached to the 
window, and are defined as the infinite boundaries of the 
finite spin chain. 

We have already introduced the idea of shrinking the 
infinite spin chain to a finite spin chain with infinite 
boundary conditions. These infinite boundary conditions 
will have to capture all the properties of the infinite sys- 
tem. Although the idea of shrinking the infinite spin 
chain is quite simple, it is more complicated to realize. 
Specifically, we need to be sure that our finite-size system 
with infinite boundary conditions will behave similarly 
to the initial infinite system. To achieve this we require 
the effective Hamiltonian representing the infinite sys- 
tem, written in the basis of the finite MPS. 



B. Effective Hamiltonian 



in the mixed canonical representation (see Fig. ^p) as 

1^-) =^...A'''-M'^'AB'*'+^B"'+^ ...|s), (2) 

s 

where A = XT and B = TX satisfy the left and right 
canonical form constraints as follows, 

^ ^..t^.. = J2 T'^^p^T"^ = I, (3) 

Si Si 

^^s.^^.t ^^p^.pLps.t ^ (4) 

Si Si 

In the above equations, and p^ are nothing but the 
right and left reduced density matrices of the spin chain 
and defined as 



X 

E 

a=l 

X 

E 



(A.)'|$^) 



(A«)'i$«; 



(5) 



(6) 



where |<i>^) and |$^) are the left and right Schmidt vec- 
tors that are orthonormal. 

The advantage of representing the MPS in the canon- 
ical form is that it not only fixes the gauge freedom in 
the MPS representation, which would otherwise cause 
numerical difficulties, but it is also very convenient for 
simplifying the computation of observables of an infinite 
system. In addition, the canonical form representation of 
IMPS is necessary in the truncation step of time evolution 
algorithms (both imaginary- and real-time evolution). 

Now, let us partition the whole infinite-size spin chain 
into three parts as illustrated in Fig.[Tj;. The middle part, 
called the window, contains N sites of the spin chain and 
the two other parts contain left and right semi-infinite 
spin chains attached to this window. Then instead of 
considering a large number of tensors outside of the win- 
dow of the iMPS we only use two matrices L and R that 



Suppose that the total Hamiltonian of the initial infi- 
nite spin chain can be decomposed into five components, 
written as 



H = Hr 



H 



LW 



(7) 



where and Hji are the Hamiltonian components for 
the left and right semi-infinite spin chain, Hi,y/ (Hwr) 
is the interaction term at the left (right) boundary of the 
window, respectively, and finally Hw is the Hamiltonian 
for the window with N sites. 

As we do not consider the whole infinite spin chain, we 
do not need the full information contained in the Hamil- 
tonian. Instead, we introduce the infinite boundary con- 
ditions to shrink the infinite chain to the finite chain. 
The Hamiltonian for this finite chain will be effectivelly 
described in the same way, as follows 



H = Ht 



H 



LW 



Hv 



H 



WR. 



Hr- 



(8) 



where the tilde symbol is added in order to distinguish 
between the effective Hamiltonian and the full Hamilto- 
nian of the system. We can see that Hw is the same 
in both Eq. [7] and Eq. [8] Our task is to find the effec- 
tive Hamiltonians of the left and right semi-infinite chain 
and their interaction components with the window (the 
components in Eq. [sjwith the tidle symbol). 

We now show the method to calculate the effective 
Hamiltonian by using spin-1 isotropic antiferromagnetic 
Heisenberg model as an example. The Hamiltonian con- 
tains nearest-neighbour interaction terms as follows 



i+li 



(9) 



where S — {S^ , S"^ , S^) is the vector containing matrices 
for the spin-1 representation of the spin algebra. 
The effective Hamiltonian can now be written as, 



H = Hl + Sl-Si + J2^^- ^'+1 +Sn-Sr + HiilO) 
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We need to find the left and right effective Hamiltonians 

Hl , Hfi and also operators Sl , Sr which are x^X matri- 
ces. The procedure to obtain the effective Hamiltonian is 
described in detail in Ref. I26|, and we now briefly review 
it here. 

Let us introduce the infinite matrix product operator 
(iMPO) which has the following form for an infinite-size 
spin chain. 



MHW) 



(11) 



where we denote |cr) = | . . . s^, Sj+i . . .) as the basis of the 
system. As the unit cell of this model contains a single 
site, the iMPO is represented by the same matrices VF*" 
repeated at every site of the chain, see Fig. [2j 




FIG. 3. (Color online): (a). Generalized transfer matrix 
for finding the left dominant eigenvector E. (b). Generalized 
transfer matrix for finding the right dominant eigenvector E. 
(c). Equation to find left dominant eigenvector, (d). Equa- 
tion to find right dominant eigenvector. 



FIG. 2. (Color online): The full Hamiltonian of the system is 
decomposed in tensor product of local matrix product opera- 
tors. 

With each type of Hamiltonian there are several ways 
to construct the iMPO, here we are using the method pro- 
posed in Ref. [2T] where all the matrices are in lower tri- 
angular forms. For the Hamiltonian described by Eq. |9j 
these matrices have the following form 



■ I 0" 

5^ 

w = sy 

S"^ 

s'^ sy s"^ I 



where I is a 3 x 3 indentity matrix. 

We now review the scheme proposed in Ref. [26 to find 
all the left effective operators; a similar scheme can be ap- 
plied for the right operators. Specifically, we need to find 
the dominant eigenvector of the transfer matrix diagra- 
matically illustrated in Fig. [3^. This dominant eigenvec- 
tor contains five components, E = {Ei, E2, E3, E4, E^). 
As we will see later, this dominant eigenvector contains 
the information of the left effective Hamiltonian that we 
need, or in DMRG terminology, E is the vector of block 
operators describing the effective Hamiltonian. However, 
as the transfer matrix is not diagonalizable, we need to 
find all the elements of E independently by employing 
the recursion relation, see Fig. [Sj:, which reads 

E^{n+1) = Tw^AEUn)) + J2 Tw,AEp{n)), (12) 

fi>a 



where we have defined 

Tx{E) = Y,{s\X\s')A''^EA', 



which is the generalized transfer operator to include a 
local operator X acting on the physical degree of freedom 
of the MPS. Note that the second term on the right-hand 



side of Eq. 12 contains only the components such that 
/3 > a, since the other terms with (3 < a are equal to 
zero as is a lower triangular matrix. Hence we can 



solve immediately the recursion relation Eq. 12 for the 
i?5 component, 

E^{n + l) = Tw.AE^{n)) 

= Y.{s\l\s')A^'^E^{n)A^. (14) 

ss' 

If the IMPS is in the canonical form, then we have 

Er,{n + l)^E^{n)^ ...^E^{Q)^l, (15) 

where I is a x x ;^ identity matrix. Moving on to £^4, we 
have 

E^in + 1) = rw«(£;4N) + 1] (sills') A^'t E^{n) A' 
^Y.{s\S^\s')A^'^A^^~Sl, (16) 

ss' 

and here the fact that the diagonal matrix element W44 = 
implies that the solution for E^ is simply a function of 
i?5 and local operators. Similarly, 

i?3(n + 1) = TwM{n))+Y,{s\Sy\s')A''^ E,{n)A'' 



Y,{s\Sy\s')A''^A' 



(17) 



(13) 



E2{n + 1) - Tw,, {E2{n)) + ^ {s\S'-'\s')A''^ E^{n) A' 
= Y.{s\S-\s')A^'^A^ ^~Sl, (18) 

ss' 

and finally, the most complicated term that contains the 
effective Hamiltonian of the left semi-infinite spin chain 
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is determined as 

Eiin + 1) = Tw,, {Ei{n)) + ^ Tw,, {Epin)) 

/3>1 

ss' 

ss' ^ y ^ 

SI 

ss' 

SI 

+ J2{s\S''\s')A''^ E^in) A^ 



(19) 



We can also write this equation in a compact form as 

Ei{n + 1)=Y, {s\lW)A''^Ei{n)A' + C, (20) 

ss' 

where C is a constant that is defined as the summation 
of last three terms in Eq. [T9j Our task is to solve Eq. [20] 
To see how this is done, let us assume the initial solution 
Ei{0) = 0. This is an arbitrary choice that has no effect 
on the final solution, up to an irrelevant constant. Then, 



Eiil) = C 

E,i2) =Tj{C)+C 

Ei{3) - Ti{Ti{C) + C) + C 

= Ti(ri(C))+T,(C) + C 



(21) 



Ei{n + 1) = TiiEiin)) + C (22) 
In general we can write the solution as follows. 



E,in)^J2^,\C) 



k=0 



= C + Tj{C) + Ti{Ti{C)) + Ti(ri(ri(C))) + . . . 

(23) 

This is the summation of a geometric series, which has 
the solution 



n-l 
fe=0 



1 - X 



In our case 



ji-i 



fc=0 



{l-T){C) ■ 



(24) 



(25) 



Notice that the spectrum of transfer matrix Tj will con- 
tain the identity I as an eigenvector with eigenvalue 1. 
Therefore, this summation will be diverging. To avoid 
this, let us decompose the summation into two terms as 



where -ffo contains all the terms that are perpendicular 
to the identity and is actually the effective Hamiltonian 
of the left semi-infinite chain, eg is a constant equal to 
the energy per site of the infinite chain. Now we can 



check the recursion relation by substituting Eq. 26 into 
Eq. [22j we have 

iJo + eo(n + 1)1 = Ti(ffo) + Ti(eoni) + C (27) 

This simplifies to a linear equation for Hq, 

{l-Ti){Ho) = C-eol (28) 

where = Tr{pC) {p is density matrix). By solving this 
linear equation we find the effective Hamiltonian and 
this completes the vector of block operators Ea {n) . Note 
that the energy per site contribution eo?iI is a constant 
shift in the energy and is therefore irrelevant for most 
purposes. 

In summary, we have explained in this section how to 
obtain the effective Hamiltonian on the left of the win- 
dow. Specifically, we have obtained — Ei and also 

operators Sl = {SI, Si, Si) = {E2,E:i,Ei). For the 
right effective Hamiltonian, a completely similar proce- 
dure is performed. In the next section we will use this 
calculation to investigate the problem of real-time evolu- 
tion of IMPS in the presence of local perturbation. 



III. APPLICATION: REAL-TIME EVOLUTION 
OF IMPS IN THE PRESENCE OF LOCAL 
PERTURBATION 

We now apply the procedure for finding the effective 
Hamiltonian and the infinite boundary conditions pro- 
posed above to study dynamical properties of an infinite 
spin chain in the presence of a local perturbation. As 
an infinite MPS will be effectively represented by a fi- 
nite MPS, we can apply a standard MPS time-evolution 
technique to study the reaction of the infinte system to a 
local perturbation. The MPS technique that we use here 
is the TEBD algorithm. 



A. Local Perturbation 

We wish to take an infinite spin chain which is in its 
ground state, and perturb locally one site. Suppose that 
we have already found the ground state of the system 
(for example, by iDMRG or iTEBD), represented by a 
translationally invariant IMPS with a one- or two-site 
unit cell; the wave function |5'gs) is written as in Eq. [2] 
Then we choose one site and perturb it locally by flipping 
the spin of that site with flipping spin operators 5"+ (flip 
spin up) or S^ (flip spin down). The system is not in the 
ground state anymore, but a mixture of excited states 
and is no longer described by a translationally invariant 
IMPS. Let us flip the spin at a certain position j in the 
chain and deflne a new state as 



Ei{n) ^ Hq + eoni 



(26) 



I«') = ^+I*GS) 



(29) 
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As a result of spin flipping, a wave packet is formed cen- 
tred at the flipped spin. As an illustration using the 
spin-1 isotropic antiferromagnetic Heisenberg model, in 
Fig. |4j we plot the local magnetization of the system after 
flipping one spin in the middle of the chain. We can see 
that a wave packet is formed with the peak located in the 
middle site. The amplitude of this wave packet decreases 
when moving away from the middle point. The width 
of the wave packet depends on the correlation length of 
the system. Note however that despite the breaking of 
translational invariance at long range, only one tensor 
of the MPS is different from that of the translationally- 
invariant ground state. 




FIG. 4. The wave packet (local magnetization) is formed 
after flipping one spin in the middle of the chain. The result 
is obtained by using a two-site translationally invariant iMPS 
for the ground state with bond dimension x — 160. 



B. Real-time evolution 

Let us now study real-time evolution of an infinite spin 
chain. The initial state of the system is a locally per- 
turbed state I^PCO)). This state will evolve in time and is 
described by the solution of the Schroedinger equation 



*(0)). 



(30) 



As mentioned above, this state can be effectively rep- 
resented by a finite MPS containing -I- 2 sites where 
the perturbed site is in the middle of the chain at site 
i = N/2, with an effective Hamiltonian H describing the 
finite system. The two boundary sites are now repre- 
sented by the boundary tensors and , which are 
the usual boundary sites of a finite MPS with dimen- 
sions 1 X X and X ^ 1 respectively, except now the local 
Hilbert space is the x-dimensional effective Hilbert space 
for the left and right semi- infinite strips. In practice we 
do not actually need the L" and tensors as these are 



identity elements; Lf = 5ai and R 







but their use 



allow us to formally write the state of the system Eq. [2] 
as a finite MPS, 



1^) 



{Si} 



..A%"R^\a,s,l3), (31) 



where |s) = |si, S2, . . . sjv), ^i"^ fullfils the condition in 
Eq. [3] and all the other tensors on the right of A matrix 



satisfy the right canonical constraint in Eq.|4] Simultane- 
uosly, the A matrix is now shifted to the link between 
and A'^^ . The time evolution process is now described by 
the unitary operator e~'^*. 

With an effective finite system representing the infi- 
nite system, we can proceed with the real-time evolution 
by employing the TEBD algorithm. Before continuing, 
we will briefiy reiterate the main features of the TEBD 
algorithm. For more details, refer to the original worlP. 
In this algorithm the time evolution operator 



iHt 



IS 



decomposed as a product of M operators e~^^^^ (where 
(5t ^ 1 is the small time step and M = t/dt). In turn, 
each term e~^^^* is decomposed into products of local 
terms by using Suzuki- Trotter decomposition. Normally, 
the Hamiltonian is written as the summation of two 
terms. With the Hamiltonian just containing nearest- 
neighbour interaction terms, we can rewrite it in the fol- 
lowing form: 



H = H, 



dd 



Hp. 



(32) 



Hp 



-^even 



where Hodd = J2odd h^'''^^^ and 
Eet,en j^'"'*"^^'- Tcrms in either Hodd 
commute with each other. However, the terms in Hodd 
do not commute with ones in Heven iu general. Then 
the first order Suzuki- Trotter decomposition of the time 
evolution operator at each time step 5t is 



odd i 



(33) 



As a consequence of the non-zero commutation relation 
between the odd and even terms of the Hamiltonian, the 
Suzuki- Trotter decomposition will produce some error on 
the order of 5t^ . However, this error can be controlled 
by using a small time step 5t or by taking high-order 
decompositioiP^. 

Here, we modify slightly the TEBD algorithm to in- 
vestigate the real-time evolution of our locally perturbed 
system. Specifically, we do not need to find the inverse 
of the A matrix after acting the two-body gate on a given 
link of two sites, but instead we use two more SVDs to 
shift the A matrix by two sites to the next update link. 
This step is also important to get an optimal truncation 
which is essential for each local update. As a price of 
implementing two extra SVDs, this step may be a little 
bit costly. However, the big advantage of doing this is 
that we can avoid the inverse of A that is numerically 
unstable. 

For convenience and clarity, we write the time evo- 
lution operator with time step 5t in the first order of 
Suzuki- Trotter decomposition 



U lU LwUwU rU Rw ■ 



(34) 
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(a) 

I ^53 ^sn J 




yl'-'i-i A"- ^'-'.-i ^'-'.-i A"- 

k 

FIG. 5. (Color online): (a). Applying the operator e~*^*' 
to the effective finite MPS. (b). Update the left and 
right tensors by contracting tensors {Ul,Ulw , I , A"^} and 
{Ur,Urw , I , A"^} . (c). Update the new tensors when apply- 
ing two-body gate it on the odd or even link. Contracting all 
the tensors involved and take the SVD of that. The bondi- 
mension will increase after taking the SVD, so we need to do 
the truncation to keep new tensors A'"' , A"'^+^ and A' in the 
desired bond dimension, (d). Shifting A' to the right (or left) 
of the updated link if sweeping direction is from left to right 
(or right to left) by taking two successive SVDs. 

The update scheme of real-time evolution at each time 
step is illustrated diagramatically in Fig. [5] for a = 6 
window. However, in principle N can be any arbitrary 
finite number. Each time step includes two successive 
sweeps: one from the left to the right and vice versa. 
Note that with our specific choice of an even number 
of sites inside the window N the interaction terms on 
the left and right sides of the window are operationally 
equivalent to the even terms. 

IV. RESULTS 

By replacing the IMPS with an effective finite MPS 
containing two boundary sites we can evolve the locally 
perturbed ground state in time. Here we present the 
results computed for the spin-1 isotropic antiferromag- 
netic Heisenberg model. The initial ground state is rep- 
resented by the IMPS with bond dimension x = 160 and 



is evolved to a state with maximum truncated bond di- 
mension xc = 200. Truncation error is approximately 
equal to 10~^. After flipping the central spin, the system 
is evolved up to time t ~ 30 where time step 6t = 0.05 is 
used for the fourth-order Suzuki- Trotter decomposition. 

A. Wave packet propagating in time 

In order to understand how the wave packet is prop- 
agating in the effective finite MPS where the infinite 
boundaries are present, we compute the local magneti- 
zations in time at each site of the spin chain. The result 
is shown in the Fig. [6] For the system with window size 
N = 60, as we can see, the wave packet at the beginning 
is formed at the middle of the chain and then spreads out. 
Importantly, what we can see here is when the wavefront 
hits the infinite boundaries, there is no back reflection or 
counter-propagating effect, it passes through the bound- 
aries. This can be verified when we look at how the 
wave packet propagates outside the window in the effec- 
tive Hilbert space. Specifically, we expand the window 
after the simulation by inserting the original orthogonal 
tensors A and B at the edges. From this we calculate 
expectation values outside the original window. From 
Fig. [6] we can see that when the window is expanded to 
window size N — 200, the wavefront moves smoothly to 
the exterior, justifying our approach. 
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FIG. 6. (Color online): Wave packet propagates in time with 
window size A^ = 60 (blue dots). The window is expanded to 
A'' = 200 (black solid lines) after simulation to see how the 
wavefront propagates beyond the infinite boundaries. 

Now, for comparison, we also plot the propagation 
of the wave packet in time with different window sizes. 
These windows have sizes fixed from the beginning of the 
real-time evolution. As we can see from Fig. [7j the wave 
packets of different window sizes are coincident with each 
other in the middle region of the plot. However, is it re- 
ally true that they are all completely coincident with each 
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other? Actually, they are not. To see the difference, let us 
plot the quantity Swi6oixi,t) = \S^,{xt,t) - SlQQ{xt,t)\ 
where Xi is the relative position of the site i^^ and the 
perturbed site which is in the middle of the chain and 
W = {60, 100, 120}. The plot in Fig. [§i corresponds to 
Xi = 0, the perturbed position, and Fig. |8]d is of the 
boundary point Xi = 30. We can see that in both plots 
the differences Swi6o{xi,t) are very small at the begin- 
ning of time evolution and then increase with time. When 
the system evolves up to t = 20, the differences increase 
approximately to 10~^ and tend to saturate in later time. 
This can be easily understood, as the wave packet has 
propagated through the infinite boundaries. 
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FIG. 7. (Color online): Wave packet propagates in time with 
(liferent window sizes which are fixed at the beginning of real- 
time evolution. The lines are distinguishable in the centre of 
the plot. 



B. Unequal time two-point correlator and spectral 
function 

Let us define an unequal time two-point correlator as 

Aix,t) = {^\S-it)S+jOM, (35) 

where the subscripts in spin-flip operators indicate po- 
sitions of the chain and xm is the middle position, |0) 
is the initial state of the system that we want to evolve. 
This equation is equivalent with 



A{x,t) 



\s-mm), 



(36) 



m which we have already 

e'"'S^{0)e-'"'^ and |V'(i)> = e-'"'S+{0)\(j)). We 
also have a phase factor appear in Eq. [36] due to \<j)) 
being the eigenvector of the Hamiltonian H correspond- 
ing to the eigenvalue Eq. Obviously, the unequal time 
two-point correlator A{x,t) can be calculated easily, as 
the time-evolved state \ip{t)) can be obtained quickly 



replaced 
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FIG. 8. (Color online): Comparisons of the differences of 
local magnetization in time at a fixed site Xi between window 
sizes N = {60, 100, 120} and window size iV = 160. (a). 
Perturbed point of chain Xi = 0; (b) Boundary point of the 
chain Xi = 30. 



from the scheme proposed above for evolving the locally 
perturbed state. 

From the unequal time two-point correlator we con- 
struct the Green's function that is defined as 



G{x,t) = -iA{x,t). 



(37) 



Fig. |9] shows the plots of the real and imaginary parts 
of the Green's function for the system with window size 
N = 60. As we can see, there are wavefronts propagating 
from the middle point toward the infinite boundaries. 
Again, there is no back reflection of the wavefront at the 
boundaries. 

By Fourier transforming of the Green's function into 
momentum and frequency spaces, we can extract the 
spectrum of the lattice system. Specifically, the Fourier 
transform of G{x,t) is 



G((7,c.) 



dt e'"* y e-*«"G(x,t). 



(38) 



For the case of the spin-1 isotropic antiferromagnetic 
Heisenberg model, the Green's function is even in x and 
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FIG. 9. (Color online): Plots of the real and imaginary parts 
of the Green's function versus time and spin chain space. 



FIG. 10. (Color online): (a). Spectral function versus mo- 
mentum and frequency for the spin-1 isotropic Heisenberg 
model, the window is = 60. (b). Spectrum viewed from 
the top when it is projected on the (w, g) plane, (c). The 
dispersion relation is derived from the maximum of the spec- 
trum. 



t, and we can simplify Eq. [38] as follows 

/oo 
dt cos Lot''^ COS qxG{x,t). (39) 

G{x, t) is a continuous function in time t. However, in 
our simulation, we have already discretized the time into 
the small time steps St. Therefore, Eq. [39]can be written 
as 

G{q,Lj) « 2 ^ cos Ljt'^ cos qxG{x,t). (40) 



The spectral function is now defined as 

S{q,U}) = --ImG{q,uj). (41) 

Note that as we have already introduced the infinite 
boundaries for our finite MPS, the wavefront can now 
propagate freely through these boundaries without any 
back refiection. In principle, in obtaining the spectral 
function we do not need to have any cutoff in time to 



keep the available data. In Fig. 10 1, we plot the spectral 
function versus momentum and frequency. In order to 
get a smooth spectrum, we have multiplied G{x,t) with 
a gaussian window function of the form exp[— 4(t/T)^] as 
introduced in Ref. [8l By viewing from the top of this 
figure, we can see the dispersion relation appears very 
clearly in Fig. 10 d. Collecting the data pairs (i^,?) that 



correspond to the maximum of spectrum and plotting 



them, we can see the dispersion relation of the system 
appears nicely in Fig. 10 The value of the gap at <? = tt 
measured in our simulation is A = 0.4105, consistent 
with the value found in Ref. [Sj 



V. CONCLUSIONS 

We have introduced the infinite boundary condition as 
a procedure for representing a finite section of a lattice 
embedded within an infinite chain. With just two bound- 
ary sites we can describe the relevant information for the 
whole semi-infinite spin chain. Therefore, instead of sim- 
ulating the iMPS, we just need the finite MPS with two 
additional effective sites. This helps to greatly reduce 
the computational cost as well as computer memory in 
simulating the infinite system where the MPS cannot be 
represented by translationally invariant tensors. After 
finding the effective Hamiltonian and operators associ- 
ated with the infinite boundary, the numerical algorithms 
we use are straightforward, variants of the well-known 
MPS/DMRG algorithms for finite size systems. Hence 
the general procedure is applicable to a wide variety of 
problems. 

As an example for possible application we considered 
the real time evolution of the ID spin-1 isotropic Heisen- 
berg model. The initial state of the system is the ground 
state where one central site is locally perturbed. As a 



10 



result, a wave packet is formed and spreads out from the 
center in time. As we have already attached the infi- 
nite boundaries to the finite system, we do not need to 
end the simulation when the wavefront hits the bound- 
aries, as useful information can be still be obtained, at 
least for short time intervals, as the degrees of freedom 
propogate into the boundary tensor. The resulting spec- 
tral function and dispersion relation compare well with 
previous investigations. The gap value we obtained com- 
pares well with that obtained in Ref. HI although smaller 
window size and longer evolution time is used. In fact, 
there is no restriction in our method which says that the 
window size must remain constant throughout the calcu- 
lation. Expanding the window size is straightforward, as 
the tensors representing the system outside the window 
are translationally invariant anyway, so additional ten- 
sors simply need to be orthogonalized and incorporated 
into the finite window. Similarly, reducing the size of the 
window is achieved by incorporating tensors from the fi- 
nite window into the infinite boundary tensors, which is 



a simple tensor contraction for the new effective Hamil- 
tonian and associated operators. The use of these tech- 
niques is described in Ref. EU 

We have described the procedure for IBC for a one- 
dimensional matrix product state, however the gen- 
eral procedure is applicable to any regular tensor net- 
work. In particular, this method is directly applicable 
to iPEPS^^^^®. This may be an effective way to obtain 
the spectral function of a 2D system, among many other 
possible applications. 
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